# R function that calculates ACPs, CACPs, differences in ACPs and CACPs,
# and direct pairwise selection probabilities for paired forced-choice
# conjoint experiments
# v8
# Flavien Ganter
# Created on July 28, 2019; last modified on February 11, 2021

# This file creates seven functions:
## 1. conjacp - A general function for calculating ACPs, DACPs, and direct 
##    pairwise preferences, that uses functions 2 and 3.
## 2. conjacp.prepdata - A function that prepares the working data set for
##    conjapc.estimation as well as all objects that conjapc.estimation requires
##    for estimating ACPs, CACPs, and direct pairwise preferences.
## 3. conjapc.estimation - A function that estimates ACPs, DACPs, and direct
##    pairwise preferences from an object generated by conjacp.prepdata.
## 4. print.conjacp - A function that specifies how conjacp objects should
##    be printed.
## 5. summary.conjacp - A function that stores the objects needed to summarize
##    conjacp objects.
## 6. print.summary.conjacp - A function that specifies how the summary of
##    conjacp objects should be printed.
## 7. conjacp.diff - A function that performs comparisons between pairs of
##    quantities of interests from conjapc.estimation.





### CONJACP ###########################################################


conjacp <- function(formula,
                    data,         
                    tasks,
                    estimand      = "acp",
                    subgroups     = NULL,
                    by            = NULL,
                    adjust        = FALSE,
                    id            = NULL,
                    level_weights = NULL,
                    condition     = NULL,
                    subset        = NULL
) {
  
  
  ### Prepare data for estimation
  if (length(by) == 1 & is.null(subgroups)) {
    subgroups <- by
    warning(paste0(by, " has been set as a subgroup variable."),
            .immediate = TRUE)
  }
  data_for_estimation <- conjacp.prepdata(formula       = formula,
                                          data          = data,         
                                          tasks         = tasks,
                                          subgroups     = subgroups,
                                          id            = id,
                                          level_weights = level_weights)
  
  
  ### Estimate quantities of interest
  out <- conjacp.estimation(conjacp.prepdata.object  = data_for_estimation,
                            estimand                 = estimand,
                            by                       = by,
                            adjust                   = adjust,
                            condition                = condition,
                            subset                   = subset)

  
  ### Output
  return(out)
  
  
}





### CONJACP.PREPDATA ###########################################################


conjacp.prepdata <- function(formula,
                             data,         
                             tasks,
                             subgroups     = NULL,
                             id            = NULL,
                             level_weights = NULL
) {
  
  
  ### PACKAGES
  
  require(sandwich)
  require(stringr)
  "%nin%" <- Negate("%in%")
  
  
  
  ### PRELIMINARY CHECKS

  
  # Outcome
  var_y <- all.vars(stats::update(formula, . ~ 0))
  if (var_y == ".") {
    stop("No outcome variable specified.")
  }
  if (class(data[, match(var_y, names(data))]) %nin% c("numeric", "integer")) {
    stop("The outcome variable should be either numeric or integer.")
  }
  if (all(names(table(data[, match(var_y, names(data))])) != c("0", "1"))) {
    stop("The outcome variable should be a dummy variable with values 0 and 1.")
  }
  outcome <- as.numeric(as.character(data[, match(var_y, names(data))]))
  
  # Attributes
  attr_x <- all.vars(stats::update(formula, 0 ~ . ))
  names(data)[which(names(data) %in% attr_x)] <-
    gsub("[\\p{P}\\p{S}\\p{Z}]", "", names(data)[which(names(data) %in% attr_x)], perl = T)
  attr_x <- gsub("[\\p{P}\\p{S}\\p{Z}]", "", attr_x, perl = T)
  if (any(lapply(data[, match(attr_x, names(data))], class) %nin% c("factor", "numeric"))) {
    stop("All attributes should be either numeric or factor.")
  }
  
  # Subgroups
  subgroups_id <- NULL
  if (length(subgroups) > 0) {
    subgroups_id        <- as.data.frame(data[, subgroups])
    names(subgroups_id) <- subgroups
    for (subgroup_temp in subgroups) {
      if (class(subgroups_id[, subgroup_temp]) != "factor") {
        stop("Subgroup variables should be factor variabes.")
      }
      if (nlevels(subgroups_id[, subgroup_temp]) != 2) {
        stop("Subgroup variables should have two, and only two, subgroups.")
      }
    }
    if (sum(is.na(subgroups_id)) != 0) {
      warning("There are missing values in the subgroup variables. Corresponding observations will be omitted in subgroup calculations.",
              immediate. = TRUE)
    }
  } else {
    subgroups_id <- data.frame(subgroup = rep(1, nrow(data)))
    subgroups    <- "subgroup"
  }
  
  # Respondents' ID
  if (is.null(id)) {
    clust <- as.vector(data[, match(tasks, names(data))])
  } else {
    clust <- as.vector(data[, match(id, names(data))])
  }
  
  # Level weights
  if (!is.null(level_weights)) {
    attr_levelweights <- c()
    for (wgt in 1:length(level_weights)) {
      if (!all(level_weights[[wgt]][[2]] >= 0)) {
        stop("Not all level weights are positive.")
      }
      level_weights[[wgt]][[1]] <-
        gsub("[\\p{P}\\p{S}\\p{Z}]", "", level_weights[[wgt]][[1]], perl = T)
      if (level_weights[[wgt]][[1]] %nin% attr_x) {
        stop(paste("Level weights are associated with unspecified attributes:", level_weights[[wgt]][[1]]))
      }
      if (class(data[[level_weights[[wgt]][[1]]]]) != "factor") {
        stop(paste("Level weights are associated with continuous attributes:", level_weights[[wgt]][[1]]))
      }
      if (nlevels(data[[level_weights[[wgt]][[1]]]]) != length(level_weights[[wgt]][[2]])) {
        stop(paste("The number of level weights do not match the number of levels:", level_weights[[wgt]][[1]]))
      }
      attr_levelweights <- c(attr_levelweights, level_weights[[wgt]][[1]])
    }
    for (attr_temp in attr_x) {
      if (class(data[[attr_temp]]) == "factor" &
          attr_temp %nin% attr_levelweights) {
        warning(paste("Level weights are only specified for some of the categorical attributes:", attr_temp))
      }
    }
  } else {
    attr_levelweights <- c()
    for (attr_temp in attr_x) {
      if (class(data[[attr_temp]]) == "factor") {
        attr_levelweights <- c(attr_levelweights, attr_temp)
        level_weights     <- list(level_weights,
                                  list(attr_temp, rep(1, nlevels(data[[attr_temp]]))))
      }
    }
    if (length(attr_levelweights) > 0) {
      level_weights <- vector(mode = "list", length = length(attr_levelweights))
      for (attr_nb in 1:length(attr_levelweights)) {
        attr_temp                <- attr_levelweights[attr_nb]
        level_weights[[attr_nb]] <- list(attr_temp, rep(1, nlevels(data[[attr_temp]])))
      }
    }
  }
  
  
  
  ### WORKING DATA SET
  
  
  ## Create base data sets
  
  # Long base data set
  data_long <- data.frame(outcome,
                          as.vector(data[, match(tasks, names(data))]),
                          clust,
                          subgroups_id,
                          data[, match(attr_x, names(data))])
  names(data_long) <- c("outcome",
                        "tasks",
                        "clust",
                        subgroups,
                        names(data)[match(attr_x, names(data))])
  
  # Store attribute labels and replace them by numbers
  attr_levels     <- vector(mode = "list", length = length(attr_x))
  attr_continuous <- NULL
  for (i in 1:length(attr_x)) {
    attr_temp <- attr_x[i]
    if (class(data_long[[attr_temp]]) == "factor") {
      attr_levels[[i]]               <- levels(data_long[[attr_temp]])
      levels(data_long[[attr_temp]]) <- 1:nlevels(data_long[[attr_temp]])
    } else {
      attr_levels[[i]] <- NA
      attr_continuous  <- c(attr_continuous, attr_temp)
    }
  }
  
  # Reshape long base data set, to wide
  data_long$option <- ave(as.numeric(data_long[["tasks"]]),
                          data_long[["tasks"]],
                          FUN = seq_along)
  data_wide <- data_long[data_long$option == 1,
                         -which(names(data_long) == "option")]
  names(data_wide)[which(names(data_wide) %in% attr_x)] <-
    paste0(names(data_wide)[which(names(data_wide) %in% attr_x)], ".1")
  data_wide_compl <- data_long[data_long$option == 2,
                               -which(names(data_long) %in% c("outcome", "option"))]
  names(data_wide_compl)[which(names(data_wide_compl) %in% attr_x)] <-
    paste0(names(data_wide_compl)[which(names(data_wide_compl) %in% attr_x)], ".2")
  data_wide <- base::merge(data_wide, data_wide_compl, by = c("tasks", "clust", subgroups))
  if (sum(complete.cases(data_wide[, -which(names(data_wide) %in% subgroups)])) <
      nrow(data_wide[, -which(names(data_wide) %in% subgroups)])) {
    warning("Some observations have been omitted due to missing values.")
    data_wide <- data_wide[complete.cases(data_wide[, -which(names(data_wide) %in% subgroups)]),]
  }
  
  
  ## Detect randomization restrictionss
  
  # Detect restricted attributes
  attr_restricted <- NULL
  formula_terms   <- terms(formula)
  if (any(attr(formula_terms, "order") >= 2)) {
    restrictions     <- colnames(attr(formula_terms, "factors")[, attr(formula_terms, "order") == 2, drop = FALSE])
    restrictions_alt <- lapply(restrictions, function(tmp) { as.formula(paste0("~", tmp)) })
    attr_restricted  <- unlist(lapply(restrictions_alt, all.vars))
    if (any(duplicated(attr_restricted))) {
      stop("All variables in constraints must be unique and constraints may only be two-way.")
    }
    attr_restricted <- unique(attr_restricted)
  }
  attr_unrestricted <- attr_x[attr_x %nin% attr_restricted]
  
  # Detect restricted levels
  levels_restricted <- NULL
  if (length(attr_restricted)) {
    levels_restricted <- vector(mode = "list", length = length(attr_restricted))
    for (attribute1 in attr_restricted) {
      
      # Identify the other attribute involved in the restriction
      restricted_pair_nb   <- which(str_detect(restrictions, attribute1))
      restricted_pair_temp <- all.vars(as.formula(paste0("~", restrictions[restricted_pair_nb])))
      attribute2           <- restricted_pair_temp[restricted_pair_temp != attribute1]
      
      # Identify restrictions
      attr_combinations               <- data.frame(data_long[[attribute1]], data_long[[attribute2]])
      attr_combinations               <- attr_combinations[!duplicated(attr_combinations),]
      attr_combinations_matrix        <- as.matrix(table(attr_combinations[,1], attr_combinations[,2]))
      indices_restrictions_attributes <- list(which(rowSums(attr_combinations_matrix) ==
                                                      min(rowSums(attr_combinations_matrix))),
                                              which(colSums(attr_combinations_matrix) ==
                                                      min(colSums(attr_combinations_matrix))))
      restrictions_attribute1         <- names(rowSums(attr_combinations_matrix)[indices_restrictions_attributes[[1]]])
      restrictions_attribute2         <- names(colSums(attr_combinations_matrix)[indices_restrictions_attributes[[2]]])
      
      # Fill list
      levels_restricted[[which(attr_restricted %in% attribute1)]] <- 
        list(attribute1 = attribute1,
             attribute2 = attribute2,
             restrictions_attribute1 = restrictions_attribute1,
             restrictions_attribute2 = restrictions_attribute2)
    }
  }
  
  
  ## Create working data set
  
  data_prep <- data.frame(outcome = data_wide[["outcome"]])
  
  # Create contrasting variables
  for (attr_temp in attr_x) {
    attribute_profile1 <- paste(attr_temp, "1", sep = ".")
    attribute_profile2 <- paste(attr_temp, "2", sep = ".")
    
    # Continuous attributes
    # Difference between profile 1 and profile 2
    if (class(data_long[[attr_temp]]) == "numeric") {
      var_temp  <- data_wide[[attribute_profile1]] - data_wide[[attribute_profile2]]
      data_prep <- cbind(data_prep, var_temp)
      colnames(data_prep)[match("var_temp", names(data_prep))] <- paste0(attr_temp, ".")
      
      # Categorical attributes
      # V_ij variables:
      ## 1 = profile 1 has the focal level and profile 2 a different level,
      ## -1 = profile 2 has the focal level and profile 1 a different level,
      ## 0 = both or none have the focal level
    } else {
      level_combinations <- combn(levels(as.factor(data_wide[[paste0(attr_temp, ".1")]])), 2)
      for (c in 1:ncol(level_combinations)) {
        var_temp  <- ifelse(as.character(data_wide[[attribute_profile1]]) == level_combinations[1,c] &
                              as.character(data_wide[[attribute_profile2]]) == level_combinations[2,c], 1,
                            ifelse(as.character(data_wide[[attribute_profile1]]) == level_combinations[2,c] &
                                     as.character(data_wide[[attribute_profile2]]) == level_combinations[1,c], -1, 0))
        data_prep <- cbind(data_prep, var_temp)
        colnames(data_prep)[match("var_temp", names(data_prep))] <-
          paste0(attr_temp, ".", level_combinations[1,c], "-", level_combinations[2,c])
      }
    }
  }
  
  # Add subgroup, cluster id, and pair id
  data_prep <- cbind(data_prep,
                     data_wide[, c(subgroups, "clust")],
                     id = 1:nrow(data_prep))
  data_prep <- data_prep[, apply(data_prep, 2, function(x) !all(x == 0))]
  
  
  
  ### OUTPUT
  
  data_for_estimation <- list(data_prep         = data_prep,
                              data_wide         = data_wide,
                              attr_x            = attr_x,
                              level_weights     = level_weights,
                              attr_levels       = attr_levels,
                              subgroups         = subgroups,
                              attr_restricted   = attr_restricted,
                              attr_unrestricted = attr_unrestricted,
                              attr_continuous   = attr_continuous,
                              levels_restricted = levels_restricted)
  return(data_for_estimation)
  
}





### CONJACP.ESTIMATION ###########################################################


conjacp.estimation <- function(conjacp.prepdata.object,
                               estimand  = "acp",
                               by        = NULL,
                               adjust    = FALSE,
                               subset    = NULL,
                               condition = NULL
) {
  
  
  ### PACKAGES
  
  require(sandwich)
  require(stringr)
  "%nin%" <- Negate("%in%")
  
  
  
  ### TRANSITION FROM CONJACP.PREPDATA
  
  
  ## Get data from conjacp.prepdata.objs
  
  data_prep         <- conjacp.prepdata.object[[1]]
  data_wide         <- conjacp.prepdata.object[[2]]
  attr_x            <- conjacp.prepdata.object[[3]]
  level_weights     <- conjacp.prepdata.object[[4]]
  attr_levels       <- conjacp.prepdata.object[[5]]
  subgroups         <- conjacp.prepdata.object[[6]]
  attr_restricted   <- conjacp.prepdata.object[[7]]
  attr_unrestricted <- conjacp.prepdata.object[[8]]
  attr_continuous   <- conjacp.prepdata.object[[9]]
  levels_restricted <- conjacp.prepdata.object[[10]]
  
  
  
  ### Sanity checks
  
  if (estimand %nin% c("acp", "dacp", "p")) {
    stop("The estimand is not recognized.")
  }
  if (!is.null(subset)) {
    if (length(subset) != 2 | !is.list(subset)) {
      stop("'subset' should be a list of two elements: the subgroup variable name and the level of that variable that should be included.")
    }
    if (subgroups %nin% subset[[1]]) {
      stop("The subset variable is not specified as a subgroup variable.")
    }
    if (subset[[2]] %nin% levels(data_prep[[subset[[1]]]])) {
      stop(paste0(subset[[2]], " is not a level of ", subset[[1]], "."))
    }
  }
  if (estimand != "dacp" & !is.null(by)) {
    stop("Subgroup calculations are only compatible with DACPs. To calculate ACPs or direct pairwise preferences, the 'by' argument should be NULL.")
  }
  if (estimand == "dacp") {
    if (!is.null(subset)) {
      stop("The 'subset' option cannot be used to calculate DACPs.")
    }
    if (is.null(subgroups) & is.null(by)) {
      stop("No subgroup variable specified.")
    }
    if (is.null(by) & length(subgroups) > 1) {
      stop("The data contains several subgroup variables. The specific subgroup variable to be used here should be specified using the 'by' argument.")
    }
    if (length(by) > 1) {
      stop("The 'by' argument should refer to one single subgroup variable.")
    }
    if (is.null(by) & length(subgroups) == 1) {
      by <- subgroups
      warning(paste0("DACPs are calculated for the subgroup variable ", subgroups, "."), .immediate = TRUE)
    }
    if (length(by) == 1 & is.null(subgroups)) {
      stop("No subgroup variables available in the transformed data. Specify the 'subgroups' argument in the conjacp.prepdata function.")
    }
  }
  if (!is.null(condition)) {
    if (length(condition) != 2 | !is.list(condition)) {
      stop("'condition' should be a list of two elements: the condition attribute name and the level(s) of that attribute that should be included.")
    }
    if (condition[[1]] %nin% attr_x) {
      stop("The condition attribute is not specified as an attribute.")
    }
    if (any(condition[[2]] %nin% attr_levels[[which(condition[[1]] == attr_x)]])) {
      stop("Not all levels specified in the condition are levels of ", condition[[1]], ".")
    }
    condition_alias <- which(attr_levels[[which(attr_x %in% condition[[1]])]] %in% condition[[2]])
  }
  
  
  ### Keep observations that belong to the specified subset
  
  if (!is.null(subset)) {
    data_prep <- data_prep[data_prep[[subset[[1]]]] == subset[[2]] & !is.na(data_prep[[subset[[1]]]]),
                           which(names(data_prep) %nin% subset[[1]])]
    data_wide <- data_wide[data_wide[[subset[[1]]]] == subset[[2]] & !is.na(data_wide[[subset[[1]]]]),
                           which(names(data_wide) %nin% subset[[1]])]
  }
  
  
  ### Keep observations that belong to the specified condition
  
  attr_cond <- ""
  if (!is.null(condition)) {
    data_prep <- data_prep[data_wide[[paste0(condition[[1]], ".1")]] %in% condition_alias &
                             data_wide[[paste0(condition[[1]], ".2")]] %in% condition_alias,]
    data_wide <- data_wide[data_wide[[paste0(condition[[1]], ".1")]] %in% condition_alias &
                             data_wide[[paste0(condition[[1]], ".2")]] %in% condition_alias,]
    attr_cond <- condition[[1]]
  }
  
  
  ### Prepare data for subgroup analysis
  
  if (length(by) == 1) {
    
    # List of subgroup variables not included
    exluded_subgroup_var <- subgroups[which(by %nin% subgroups)]
    
    # Remove NAs, keep the focal subgroup variable, and rename levels - data_prep
    data_prep <- data_prep[!is.na(data_prep[[by]]),
                           which(names(data_prep) %nin% exluded_subgroup_var)]
    names(data_prep)[which(names(data_prep) == by)] <- "subgroup"
    subgroup_levels <- levels(data_prep$subgroup)
    levels(data_prep$subgroup) <- c("1", "2")
    
    # Remove NAs, keep the focal subgroup variable, and rename levels - data_wide
    data_wide <- data_wide[!is.na(data_wide[[by]]),
                           which(names(data_wide) %nin% exluded_subgroup_var)]
    names(data_wide)[which(names(data_wide) == by)] <- "subgroup"
    subgroup_levels <- levels(data_wide$subgroup)
    levels(data_wide$subgroup) <- c("1", "2")
    
  } else {
    if (length(subgroups) > 0) {
      data_prep <- data_prep[, which(names(data_prep) %nin% subgroups)]
    }
    data_prep$subgroup <- 1
    data_wide$subgroup <- 1
  }
  
  
  ### ESTIMATION
    
    
  ## Estimate OLS regression for each attribute
    
  # Lists to store model outputs
  pairwise_preferences <- c()
  list_vcov <- vector(mode = "list", length = (length(attr_unrestricted[attr_unrestricted != attr_cond]) +
                                                 2 * length(attr_restricted[attr_restricted != attr_cond])))
  list_vcov_indices <- 0
    
  # Loop over attributes
  for (attr_temp in attr_x[attr_x != attr_cond]) {
    
    # Loop over subgroups
    for (sub in unique(data_prep$subgroup)) {
      
      # Prepare data set
      data_model <- as.data.frame(data_prep[data_prep$subgroup == sub,])
      names(data_model)[which(str_extract(names(data_model), ".+?(?=\\.)") != attr_temp)] <- 
        paste0("c_",
               names(data_model)[which(str_extract(names(data_model), ".+?(?=\\.)") != attr_temp)],
               ".",
               attr_temp)
      if (attr_temp %in% attr_continuous) {
        names(data_model)[which(names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))] <- 
          paste0(names(data_model)[which(names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))],
                 ".",
                 sub)
      } else {
        names(data_model)[which(names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))] <- 
          paste0(names(data_model)[which(names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))],
                 "..",
                 sub)
      }
      variables_names_reorder <- c(names(data_model)[which(names(data_model) %in% c("outcome", "clust", "id", "subgroup"))],
                                   names(data_model)[which(substr(names(data_model), 1, 2) != "c_" &
                                                             names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))],
                                   names(data_model)[which(substr(names(data_model), 1, 2) == "c_" &
                                                             names(data_model) %nin% c("outcome", "clust", "id", "subgroup"))])
      data_model <- data_model[, variables_names_reorder]
      
      # Adjustement
      if (adjust == FALSE) {
        data_model <- data_model[, -which(substr(colnames(data_model), 1, 2) == "c_")]
      }
      
      # For unconstrained attributes, run a unique regression
      if (attr_temp %in% attr_unrestricted) {
        
        # Prepare data
        cluster_temp <- data_model[["clust"]]
        id_temp      <- data_model[["id"]]
        data_model   <- data_model[, -match(c("clust", "id", "subgroup"), names(data_model))]
        names(data_model)[which(names(data_model) != "outcome")] <-
          paste0(names(data_model)[which(names(data_model) != "outcome")], "-0")
        
        # Run model
        model_temp <- lm(outcome ~ ., data = data_model)
        
          # Warnings
          if (length(unique(cluster_temp)) < 42) {
            warning(paste0("Choices of less than 42 respondents are used in the estimation of the quantities of interest for ", attr_temp, ". Standard errors might be unreliable."),
                    immediate. = TRUE)
          }
          if (model_temp$rank == nrow(data_model)) {
            warning(paste0("In the estimation of the quantities of interest for ", attr_temp, ", the number of observations equals the number of estimated parameters. A slight correction of the degree of freedom adjustment is performed; standard errors might be unreliable."),
                    immediate. = TRUE)
          }
          if (nrow(data_model) < 100) {
            warning(paste0("The quantities of interest for ", attr_temp, "are estimated on less than 100 choices."),
                    immediate. = TRUE)
          }
        
        # Store coefficients
        coef_temp <- coefficients(model_temp)
        coef_temp <- coef_temp[-which(names(coef_temp) == "(Intercept)")]
        if (adjust == TRUE) {
          coef_temp <- coef_temp[-which(substr(names(coef_temp), 2, 3) == "c_")]
        }
        if (sum(is.na(coef_temp))) {
          stop("Some direct pairwise comparisons cannot be estimated from the data, either because they do not exist in the pool of profiles, or because those in the pool are collinear with other comparisons.")
        }
        names(coef_temp)     <- gsub("`", "", names(coef_temp), fixed = TRUE)
        pairwise_preferences <- c(pairwise_preferences, coef_temp)
        
        # Store elements for variance-covariance matrix
        list_vcov_indices              <- list_vcov_indices + 1
        list_vcov[[list_vcov_indices]] <- list(model_temp, id_temp, cluster_temp)
        
      } else {
      # For constrainted attributes, run two regressions:
        ## (1) one regression conditional on the unrestricted levels of the second attribute, for all levels
        ## to be fully comparable;
        ## (3) one regression conditional on comparable pairs, for unrestricted levels (for restricted levels, the CACP is
        ## the same as in (1)---à la AMCE.
        
        attribute1 <- attr_temp
        
        for (r in 1:2) {
          
          # Copy model data set
          data_model_cacp <- data_model
          
          # Construct restricted data sets for CACPs
          
            # Extract restrictions
            attribute2              <- levels_restricted[[which(attr_restricted %in% attribute1)]]$attribute2
            attribute1_restrictions <- levels_restricted[[which(attr_restricted %in% attribute1)]]$restrictions_attribute1
            attribute2_restrictions <- levels_restricted[[which(attr_restricted %in% attribute1)]]$restrictions_attribute2
            
            # Add guide variables
            data_model_cacp <- cbind(data_model_cacp,
                                     data_wide[data_wide$subgroup == sub,
                                               c(paste0(attribute1, ".1"),
                                                 paste0(attribute1, ".2"),
                                                 paste0(attribute2, ".1"),
                                                 paste0(attribute2, ".2"))])
            names(data_model_cacp)[c((length(names(data_model_cacp))-3):
                                       length(names(data_model_cacp)))] <- c("attribute1.1", "attribute1.2",
                                                                             "attribute2.1", "attribute2.2")
            
            # Remove observations
            if (r == 1) {
              data_model_cacp <- data_model_cacp[data_model_cacp$attribute2.1 %nin% attribute2_restrictions &
                                                   data_model_cacp$attribute2.2 %nin% attribute2_restrictions,]
            } else {
              data_model_cacp <- data_model_cacp[(((data_model_cacp$attribute1.1 %in% attribute1_restrictions &
                                                      data_model_cacp$attribute1.2 %nin% attribute1_restrictions) |
                                                     (data_model_cacp$attribute1.1 %nin% attribute1_restrictions &
                                                        data_model_cacp$attribute1.2 %in% attribute1_restrictions)) &
                                                    data_model_cacp$attribute2.1 %nin% attribute2_restrictions &
                                                    data_model_cacp$attribute2.2 %nin% attribute2_restrictions) |
                                                   (data_model_cacp$attribute1.1 %nin% attribute1_restrictions &
                                                      data_model_cacp$attribute1.2 %nin% attribute1_restrictions) |
                                                   (data_model_cacp$attribute1.1 %in% attribute1_restrictions &
                                                      data_model_cacp$attribute1.2 %in% attribute1_restrictions),]
            }
            data_model_cacp <- data_model_cacp[, apply(data_model_cacp, 2, function(x) !all(x == 0))]
          
            # Remove guide variables
            data_model_cacp <- data_model_cacp[, -match(c("attribute1.1", "attribute1.2", "attribute2.1", "attribute2.2"),
                                                        names(data_model_cacp))]
          
          # Store clust and id, and remove them from the data set
          cluster_temp    <- data_model_cacp[["clust"]]
          id_temp         <- data_model_cacp[["id"]]
          data_model_cacp <- data_model_cacp[, -match(c("clust", "id", "subgroup"), names(data_model_cacp))]
          names(data_model_cacp)[which(names(data_model_cacp) != "outcome")] <-
            paste0(names(data_model_cacp)[which(names(data_model_cacp) != "outcome")], "-", r)
          
          # Run model
          model_temp <- lm(outcome ~ ., data = data_model_cacp)
          
            # Warnings
            if (length(unique(cluster_temp)) < 42) {
              warning(paste0("Choices of less than 42 respondents are used in the estimation of the conditional (c", r, ") quantities of interest for ", attr_temp, ". Standard errors might be unreliable."),
                      immediate. = TRUE)
            }
            if (model_temp$rank == nrow(data_model)) {
              warning(paste0("In the estimation of the conditional (c", r, ") quantities of interest for ", attr_temp, ", the number of observations equals the number of estimated parameters. A slight correction of the degree of freedom adjustment is performed; standard errors might be unreliable."),
                      immediate. = TRUE)
            }
            if (nrow(data_model) < 100) {
              warning(paste0("The conditional (c", r, ") quantities of interest for ", attr_temp, "are estimated on less than 100 choices."),
                      immediate. = TRUE)
            }
          
          # Store coefficients
          coef_temp <- coefficients(model_temp)
          coef_temp <- coef_temp[-which(names(coef_temp) == "(Intercept)")]
          if (adjust == TRUE) {
            coef_temp <- coef_temp[-which(substr(names(coef_temp), 2, 3) == "c_")]
          }
          if (sum(is.na(coef_temp))) {
            stop("Some direct pairwise comparisons cannot be estimated from the data, either because they do not exist in the pool of profiles, or because those in the pool are collinear with other comparisons.")
          }
          names(coef_temp)     <- gsub("`", "", names(coef_temp), fixed = TRUE)
          pairwise_preferences <- c(pairwise_preferences, coef_temp)
          
          # Store elements for variance-covariance matrix
          list_vcov_indices              <- list_vcov_indices + 1
          list_vcov[[list_vcov_indices]] <- list(model_temp, id_temp, cluster_temp)
          
        }
        
      }
      
    }
  }
  
  
  ## Estimation of the full variance-covariance matrix
  
  # Function: clustered variance-covariance matrix
  vcov.cluster <- function(model, cluster) {
    require(sandwich)
    
    # Residuals
    uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum))
    
    # Degree of freedom correction
    N    <- length(cluster)        
    r    <- model$rank
    N_id <- length(unique(cluster))
    if (N != r) {
      dfc <- (N_id/(N_id - 1)) * ((N - 1)/(N - r))
    } else {
      # Change the degree of freedom adjustment if N = r --- standard errors might not be reliable
      dfc <- (N_id/(N_id - 1)) * ((N + 1 - 1)/(N + 1 - r))
    }
    
    clustered_vcov <- dfc * sandwich(model, meat = crossprod(uj)/N)
    return(clustered_vcov)
  }
  
  # Function: between-model coefficient covariance
  cov.2models <- function(list_model1, list_model2) {
    require(sandwich)
    
    # Unpack list
    model1   <- list_model1[[1]]
    id1      <- list_model1[[2]]
    cluster1 <- list_model1[[3]]
    model2   <- list_model2[[1]]
    id2      <- list_model2[[2]]
    cluster2 <- list_model2[[3]]
    
    # Observations used in both models
    id      <- id1[id1 %in% id2]
    cluster <- cluster1[id1 %in% id2]
    
    # Model matrices
    X_1 <- model.matrix(model1)
    X_2 <- model.matrix(model2)
    
      # Remove collinear vectors, if any
      indices_na1 <- which(is.na(coefficients(lm(rep(1, dim(X_1)[1])~0+X_1))))
      indices_na2 <- which(is.na(coefficients(lm(rep(1, dim(X_2)[1])~0+X_2))))
      if (length(indices_na1)) {
        X_1 <- X_1[,-indices_na1]
      }
      if (length(indices_na2)) {
        X_2 <- X_2[,-indices_na2]
      }
    
    # Inverse Jacobian matrices
    XX1_1 <- solve(crossprod(X_1))#, tol = 1e-100)
    XX1_2 <- solve(crossprod(X_2))#, tol = 1e-100)
    
    # Residuals
    estfun_1 <- estfun(model1)[id1 %in% id,]
    estfun_2 <- estfun(model2)[id2 %in% id,]
    uj_1     <- apply(estfun_1, 2, function(x) tapply(x, as.character(cluster), sum))
    uj_2     <- apply(estfun_2, 2, function(x) tapply(x, as.character(cluster), sum))
    
    # Degree of freedom correction
    N    <- 2 * length(id)          
    r    <- model1$rank + model2$rank
    N_id <- length(unique(cluster))
    dfc  <- (N_id/(N_id - 1)) * ((N - 1)/(N - r))
    
    # Obtain covariances
    if (length(uj_1) == 0) {
      cov12 <- matrix(0, nrow = ncol(X_1), ncol = ncol(X_2))
    } else {
      cov12 <- dfc * XX1_1 %*% crossprod(uj_1,uj_2) %*% XX1_2
    }
    return(cov12)
    
  }
  
  # Function: get full variance-covariance matrix
  get.full.vcov <- function(models) {
    
    # Get dimension of the final vcov matrix
    rank <- c()
    for (m in 1:length(models)) {
      rank[m] <- models[[m]][[1]]$rank
    }
    dim_vcov <- sum(rank)
    
    # Create empty vcov matrix
    full_vcov <- matrix(NA, nrow = dim_vcov, ncol = dim_vcov)
    
    # Calculate diagonal matrices
    full_vcov_names <- c()
    for (m in 1:length(models)) {
      coord_no <- sum(rank[0:(m-1)])
      full_vcov[(coord_no+1):(coord_no+rank[m]),
                (coord_no+1):(coord_no+rank[m])] <- vcov.cluster(models[[m]][[1]], models[[m]][[3]])
      full_vcov_names <- c(full_vcov_names,
                           gsub("`", "", names(models[[m]][[1]]$coefficients[!is.na(models[[m]][[1]]$coefficients)]), fixed = TRUE))
    }
    dimnames(full_vcov) <- list(full_vcov_names, full_vcov_names)
    
    # Calculate covariance matrices
    all_pairs <- combn(length(models), 2)
    for (pair in 1:ncol(all_pairs)) {
      
      # Set indices
      i <- all_pairs[1, pair]
      j <- all_pairs[2, pair]
      
      # Position in full matrix
      coord_no_1 <- sum(rank[0:(i-1)])
      coord_no_2 <- sum(rank[0:(j-1)])
      
      # Calculate covariance
      cov_temp <- cov.2models(models[[i]], models[[j]])
      
      # Plug covariance in full matrix
      full_vcov[(coord_no_2+1):(coord_no_2+rank[j]),
                (coord_no_1+1):(coord_no_1+rank[i])] <- t(cov_temp)
      full_vcov[(coord_no_1+1):(coord_no_1+rank[i]),
                (coord_no_2+1):(coord_no_2+rank[j])] <- cov_temp
      
    }
    
    # Out
    return(full_vcov)
    
  }
  
  # Calculate full variance-covariance matrix
  full_vcov <- get.full.vcov(list_vcov)
  
  # Remove rows and cols related to intercepts and adjusting variables
  names_vcov <- dimnames(full_vcov)[[1]]
  full_vcov  <- full_vcov[-which(names_vcov == "(Intercept)" |
                                   substr(names_vcov, 1, 2) == "c_"),
                          -which(names_vcov == "(Intercept)" |
                                   substr(names_vcov, 1, 2) == "c_")]
  
  
  
  ### CALCULATE QOI
    
    
  # ACPs, CACPs, and DACPs
  if (estimand %in% c("acp", "dacp")) {
    
    ## Store coefficients in a data frame
    
    # Create empty data frame
    coef_df <- data.frame(name      = names(pairwise_preferences),
                          value     = pairwise_preferences,
                          attribute = NA,
                          level1    = NA,
                          level2    = NA,
                          subgroup  = NA,
                          cacp      = NA)
    coef_df$name <- as.character(coef_df$name)
    
    # Fill table
    coef_df$attribute <- str_extract(coef_df$name, ".+?(?=\\.)")
    coef_df$level1    <- str_extract(coef_df$name, "(?<=\\.).+?(?=\\-)")
    coef_df$level2    <- str_extract(coef_df$name, "(?<=\\-).+?(?=\\.\\.)")
    coef_df$subgroup  <- str_extract(coef_df$name, "(?<=\\.\\.).+?(?=\\-)")
    coef_df$cacp      <- substr(coef_df$name, nchar(coef_df$name), nchar(coef_df$name))
    
    
    ## Scalar calculation
    
      # ACPs, CACPs, and DACPs are linear combinations of the direct pairwise preferences previously
      # estimated. This section constructs the coefficients of the linear combinations, which are
      # all be stored in a matrix. The quantities of interest is obtained by multiplying the
      # coefficient matrix and the vector of pairwise preferences. The variance-covariance matrix
      # is claculated by multiplying the coefficient matrix and the full variance-covariance matrix
      # of the direct pairwise preferences.
    
    # Create empty object to store scalars
    scalars <- NULL
    
    # Loop over attributes
    for (attr_temp in attr_x[attr_x != attr_cond]) {
      
      # Continuous attributes
      if (attr_temp %in% attr_continuous) {
        
        # Determine scalars
        scalars_temp <- ifelse(coef_df[["attribute"]] == attr_temp &
                                 coef_df[["subgroup"]] == 1, 1, 0)
        
        # Adjust if there are several subgroups
        if (estimand == "dacp") {
          scalars_temp_2 <- ifelse(coef_df[["attribute"]] == attr_temp &
                                     coef_df[["subgroup"]] == 2, 1, 0)
          scalars_temp   <- scalars_temp - scalars_temp_2
        }
        
        # Store the vector of scalars in the full matrix
        scalars <- as.data.frame(cbind(scalars, scalars_temp))
        colnames(scalars)[match("scalars_temp", names(scalars))] <- attr_temp
        
      # Categorical, unrestricted attributes
      } else if (attr_temp %in% attr_unrestricted) {
        
        # Extract level weights
        level_weights_temp <- level_weights[[which(lapply(level_weights, function(l) l[[1]]) %in% attr_temp)]][[2]]
        
        # Determine variable levels
        attr_temp_levels <- as.character(1:length(level_weights_temp))
        
        
        # Loop over attribute levels
        for (level_temp in attr_temp_levels) {
          
          # Determine the non-normalized value of each scalar
          scalars_temp <- ifelse(coef_df[["attribute"]] == attr_temp &
                                   coef_df[["level1"]] == level_temp &
                                   coef_df[["cacp"]] == 0 &
                                   coef_df[["subgroup"]] == 1,
                                 level_weights_temp[as.numeric(coef_df$level2)],
                                 ifelse(coef_df[["attribute"]] == attr_temp &
                                          coef_df[["level2"]] == level_temp &
                                          coef_df[["cacp"]] == 0 &
                                          coef_df[["subgroup"]] == 1,
                                        -1 * level_weights_temp[as.numeric(coef_df$level2)],
                                        0))
          
          # Take the difference if there are several subgroups
          if (estimand == "dacp") {
            scalars_temp_2 <- ifelse(coef_df[["attribute"]] == attr_temp &
                                       coef_df[["level1"]] == level_temp &
                                       coef_df[["cacp"]] == 0 &
                                       coef_df[["subgroup"]] == 2,
                                     level_weights_temp[as.numeric(coef_df$level2)],
                                     ifelse(coef_df[["attribute"]] == attr_temp &
                                              coef_df[["level2"]] == level_temp &
                                              coef_df[["cacp"]] == 0 &
                                              coef_df[["subgroup"]] == 2,
                                            -1 * level_weights_temp[as.numeric(coef_df$level2)],
                                            0))
            scalars_temp <- scalars_temp - scalars_temp_2
          }
          
          # Normalize scalars
          sum_level_weights <- sum(level_weights_temp[-as.numeric(level_temp)])
          scalars_temp      <- 1/sum_level_weights * scalars_temp
          
          # Store the vector of scalars in the full matrix
          scalars <- as.data.frame(cbind(scalars, scalars_temp))
          colnames(scalars)[match("scalars_temp", names(scalars))] <- paste0(attr_temp, ".", level_temp)
          
        }
        
      # Categorical, restricted attributes
      } else {
        
        # Extract level weights
        level_weights_temp <- level_weights[[which(lapply(level_weights, function(l) l[[1]]) %in% attr_temp)]][[2]]
        
        # Loop over CACPs
        for (r in 1:2) {
          
          # Determine variable levels
          attr_temp_levels <- as.character(1:length(level_weights_temp))
          
          # Loop over attribute levels
          for (level_temp in attr_temp_levels) {
            
            # Determine the non-normalized value of each scalar
            scalars_temp <- ifelse(coef_df[["attribute"]] == attr_temp &
                                     coef_df[["level1"]] == level_temp &
                                     coef_df[["cacp"]] == r &
                                     coef_df[["subgroup"]] == 1,
                                   level_weights_temp[as.numeric(coef_df$level2)],
                                   ifelse(coef_df[["attribute"]] == attr_temp &
                                            coef_df[["level2"]] == level_temp &
                                            coef_df[["cacp"]] == r &
                                            coef_df[["subgroup"]] == 1,
                                          -1 * level_weights_temp[as.numeric(coef_df$level2)],
                                          0))
            
            # Take the difference if there are several subgroups
            if (estimand == "dacp") {
              scalars_temp_2 <- ifelse(coef_df[["attribute"]] == attr_temp &
                                         coef_df[["level1"]] == level_temp &
                                         coef_df[["cacp"]] == r &
                                         coef_df[["subgroup"]] == 2,
                                       level_weights_temp[as.numeric(coef_df$level2)],
                                       ifelse(coef_df[["attribute"]] == attr_temp &
                                                coef_df[["level2"]] == level_temp &
                                                coef_df[["cacp"]] == r &
                                                coef_df[["subgroup"]] == 2,
                                              -1 * level_weights_temp[as.numeric(coef_df$level2)],
                                              0))
              scalars_temp <- scalars_temp - scalars_temp_2
            }
            
            # Normalize scalar
            sum_level_weights <- sum(level_weights_temp[-as.numeric(level_temp)])
            scalars_temp      <- 1/sum_level_weights * scalars_temp
            
            # Store the vector of scalars in the full matrix
            scalars <- as.data.frame(cbind(scalars, scalars_temp))
            colnames(scalars)[match("scalars_temp", names(scalars))] <- paste0(attr_temp, ".", level_temp, "-c", r)
            
          }
          
        }
        
      }
      
    }
    
  }
  
  
  ## Get quantities of interest
  
  # Function: change levels back
  get_levels_back <- function(coef_name) {
    
    if (estimand %in% c("acp", "dacp")) {
    
      # Decompose coefficient name
      attribute_name <- str_extract(coef_name, ".+?(?=\\.)")
      
      if (is.na(attribute_name)) {
        
        new_coef_name <- coef_name
        
      } else {
      
        cacp_indicator <- NULL
        if (substr(coef_name, nchar(coef_name)-2, nchar(coef_name)) %in% c("-c1", "-c2", "-c3")) {
          cacp_indicator <- paste0("..", substr(coef_name, nchar(coef_name)-1, nchar(coef_name)))
          coef_name      <- substr(coef_name, 1, nchar(coef_name)-3)
        }
        focal_level <- as.numeric(str_extract(coef_name, "(?<=\\.).+"))
        
        # Compose new coefficient name
        attr_index    <- which(attr_x %in% attribute_name)
        level_text    <- attr_levels[[attr_index]][focal_level]
        new_coef_name <- paste0(attribute_name, ".", level_text, cacp_indicator)
      
      }
    
    } else {
      
      # Decompose coefficient name
      attribute_name <- str_extract(coef_name, ".+?(?=\\.)")
      
      if (attribute_name %in% attr_continuous) {
        
        new_coef_name <- attribute_name
        
      } else {
      
        cacp_indicator <- NULL
        if (substr(coef_name, nchar(coef_name), nchar(coef_name)) %in% c("1", "2", "3")) {
          cacp_indicator <- paste0("..c", substr(coef_name, nchar(coef_name), nchar(coef_name)))
        }
        coef_name <- substr(coef_name, 1, nchar(coef_name)-5)
        level1    <- as.numeric(str_extract(coef_name, "(?<=\\.).+?(?=\\-)"))
        level2    <- as.numeric(str_extract(coef_name, "(?<=\\-).+"))
        
        # Compose new coefficient name
        attr_index    <- which(attr_x %in% attribute_name)
        level1_text   <- attr_levels[[attr_index]][level1]
        level2_text   <- attr_levels[[attr_index]][level2]
        new_coef_name <- paste0(attribute_name, ".", level1_text, "-", level2_text, cacp_indicator)
        
      }
      
    }
    
    return(new_coef_name)
    
  }
  
  # ACPs, CACPs, and DACPs
  if (estimand %in% c("acp", "dacp")) {
  
    # Convert the scalar data frame into a matrix object
    scalars <- as.matrix(scalars)
    
    # Get coefficient names
    names_acp <- sapply(colnames(scalars), get_levels_back)
    
    # Indices for subsets of estimates - original
    if (length(attr_unrestricted) > 0) {
      indices_acp          <- which(substr(names_acp, nchar(names_acp)-3, nchar(names_acp)-1) != "..c")
    } else {
      indices_acp          <- NULL
    }
    if (length(attr_restricted) > 0) {
      indices_cacp1        <- which(substr(names_acp, nchar(names_acp)-3, nchar(names_acp)) == "..c1")
      indices_cacp2        <- which(substr(names_acp, nchar(names_acp)-3, nchar(names_acp)) == "..c2")
    } else {
      indices_cacp1        <- NULL
      indices_cacp2        <- NULL
    }
    
    # Point estimates
    estimates_all        <- t(scalars) %*% pairwise_preferences
    estimates_all        <- as.vector(estimates_all)
    estimates            <- estimates_all[c(indices_acp, indices_cacp1)]
    names(estimates)     <- c(names_acp[indices_acp], str_extract(names_acp[indices_cacp1], ".+?(?=\\.\\.)"))
    estimates_alt        <- estimates_all[c(indices_acp, indices_cacp2)]
    names(estimates_alt) <- c(names_acp[indices_acp], str_extract(names_acp[indices_cacp2], ".+?(?=\\.\\.)"))
    
    # Variance-covariance matrix
    vcov_all           <- t(scalars) %*% full_vcov %*% scalars
    vcov               <- vcov_all[c(indices_acp, indices_cacp1), c(indices_acp, indices_cacp1)]
    dimnames(vcov)     <- list(c(names_acp[indices_acp], str_extract(names_acp[indices_cacp1], ".+?(?=\\.\\.)")),
                               c(names_acp[indices_acp], str_extract(names_acp[indices_cacp1], ".+?(?=\\.\\.)")))
    vcov_alt           <- vcov_all[c(indices_acp, indices_cacp2), c(indices_acp, indices_cacp2)]
    dimnames(vcov_alt) <- list(c(names_acp[indices_acp], str_extract(names_acp[indices_cacp2], ".+?(?=\\.\\.)")),
                               c(names_acp[indices_acp], str_extract(names_acp[indices_cacp2], ".+?(?=\\.\\.)")))
    
    # Indices for subsets of estimates - out
    if (length(attr_unrestricted) > 0) {
      indices_acp         <- which(str_extract(names(estimates), ".+?(?=\\.)") %in% attr_unrestricted)
    }
    if (length(attr_restricted) > 0) {
      indices_cacp        <- which(str_extract(names(estimates), ".+?(?=\\.)") %in% attr_restricted)
      indices_cacp_byattr <- vector(mode = "list", length = length(attr_restricted))
      for (attr_temp in attr_restricted) {
        attr_nb                         <- which(attr_restricted %in% attr_temp)
        indices_cacp_byattr_temp        <- which(str_extract(names(estimates), ".+?(?=\\.)") == attr_temp)
        indices_cacp_byattr[[attr_nb]]  <- list(attribute = attr_temp,
                                                indices = indices_cacp_byattr_temp[
                                                  indices_cacp_byattr_temp %in% indices_cacp])
      }
    } else {
      indices_cacp        <- NULL
      indices_cacp_byattr <- NULL
    }
    
  # P
  } else {
    
    # Get coefficient names
    names_p <-  sapply(names(pairwise_preferences), get_levels_back)
    
    # Indices for subsets of estimates - original
    if (length(attr_unrestricted) > 0) {
      indices_acp          <- which(substr(names_p, nchar(names_p)-3, nchar(names_p)-1) != "..c")
    } else {
      indices_acp          <- NULL
    }
    if (length(attr_restricted) > 0) {
      indices_cacp1        <- which(substr(names_p, nchar(names_p)-3, nchar(names_p)) == "..c1")
      indices_cacp2        <- which(substr(names_p, nchar(names_p)-3, nchar(names_p)) == "..c2")
    } else {
      indices_cacp1        <- NULL
      indices_cacp2        <- NULL
    }
    
    # Point estimates
    estimates_all        <- pairwise_preferences
    estimates            <- estimates_all[c(indices_acp, indices_cacp1)]
    names(estimates)     <- c(names_p[indices_acp], str_extract(names_p[indices_cacp1], ".+?(?=\\.\\.)"))
    estimates_alt        <- estimates_all[c(indices_acp, indices_cacp2)]
    names(estimates_alt) <- c(names_p[indices_acp], str_extract(names_p[indices_cacp2], ".+?(?=\\.\\.)"))
    
    # Variance-covariance matrix
    vcov_all           <- full_vcov
    vcov               <- vcov_all[c(indices_acp, indices_cacp1), c(indices_acp, indices_cacp1)]
    dimnames(vcov)     <- list(c(names_p[indices_acp], str_extract(names_p[indices_cacp1], ".+?(?=\\.\\.)")),
                               c(names_p[indices_acp], str_extract(names_p[indices_cacp1], ".+?(?=\\.\\.)")))
    vcov_alt           <- vcov_all[c(indices_acp, indices_cacp2), c(indices_acp, indices_cacp2)]
    dimnames(vcov_alt) <- list(c(names_p[indices_acp], str_extract(names_p[indices_cacp2], ".+?(?=\\.\\.)")),
                               c(names_p[indices_acp], str_extract(names_p[indices_cacp2], ".+?(?=\\.\\.)")))
    
    # Indices for subsets of estimates - out
    if (length(attr_unrestricted) > 0) {
      indices_acp         <- which(str_extract(names(estimates), ".+?(?=\\.)") %in% attr_unrestricted)
    }
    if (length(attr_restricted) > 0) {
      indices_cacp        <- which(str_extract(names(estimates), ".+?(?=\\.)") %in% attr_restricted)
      indices_cacp_byattr <- vector(mode = "list", length = length(attr_restricted))
      for (attr_temp in attr_restricted) {
        attr_nb                         <- which(attr_restricted %in% attr_temp)
        indices_cacp_byattr_temp        <- which(str_extract(names(estimates), ".+?(?=\\.)") == attr_temp)
        indices_cacp_byattr[[attr_nb]]  <- list(attribute = attr_temp,
                                                indices = indices_cacp_byattr_temp[
                                                  indices_cacp_byattr_temp %in% indices_cacp])
      }
    } else {
      indices_cacp        <- NULL
      indices_cacp_byattr <- NULL
    }
    
  }
    
    
  
  ### OUTPUT
    
  # Subgroup parameters
  subgroups_parameters <- NULL
  if (estimand == "dacp") {
    subgroups_parameters <- list(var = by,
                                 ref = subgroup_levels[1],
                                 oth = subgroup_levels[2])
  }
  
  # Condition parameters
  if (attr_cond == "") { attr_cond <- NULL }
  
  # Out
  out <- list(estimates                   = estimates,
              vcov                        = vcov,
              estimates_alt               = estimates_alt,
              vcov_alt                    = vcov_alt,
              estimand                    = estimand,
              subgroups_parameters        = subgroups_parameters,
              n_profiles                  = nrow(data_prep) * 2,
              n_tasks                     = nrow(data_prep),
              n_respondents               = length(unique(data_prep$clust)),
              subset                      = subset,
              attr_list                   = attr_x,
              attr_restricted             = attr_restricted,
              levels_restricted           = levels_restricted,
              attr_levels                 = attr_levels,
              condition                   = condition,
              indices                     = list(acp = indices_acp,
                                                 cacp = indices_cacp,
                                                 cacp_byattr = indices_cacp_byattr))

  class(out) <- "conjacp"
  return(out) 
  
}





### PRINT ###########################################################


print.conjacp <- function(x,
                          cacp = "conditional",
                          digits = max(3L, getOption("digits") - 3L)
) {
  if (length(x$indices$acp)) {
    if (x$estimand == "acp") {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nACPs:\n")
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nACPs: (unrestricted attributes)\n")
      }
    } else if (x$estimand == "dacp") {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nDACPs,",
            paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":\n"))
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nDACPs,",
            paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":  (unrestricted attributes)\n"))
      }
    } else {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nDirect Pairwise Preferences:\n")
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nDirect Pairwise Preferences: (unrestricted attributes)\n")
      }
    }
    print.default(format(x$estimates[x$indices$acp], digits = digits), print.gap = 2L, quote = FALSE)
  }
  if (length(x$indices$cacp_byattr)) {
    if (x$estimand == "acp") {
      cat("\nCACPs:\n")
    } else if (x$estimand == "dacp") {
      cat("\nDCACPs,",
          paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":\n"))
    } else {
      cat("\nDirect Pairwise Preferences:\n")
    }
    if (cacp == "conditional") {
      for (attr_nb in 1:length(x$indices$cacp_byattr)) {
        attribute1              <- x$indices$cacp_byattr[[attr_nb]][[1]]
        attribute2              <- x$levels_restricted[[which(x$attr_restricted %in% attribute1)]]$attribute2
        attribute1_restrictions <- as.numeric(x$levels_restricted[[which(x$attr_restricted %in% attribute1)]]$restrictions_attribute2)
        conditions_cacp1        <- x$attr_levels[[which(x$attr_list %in% attribute2)]][-attribute1_restrictions]
        conditions_cacp2        <- x$attr_levels[[which(x$attr_list %in% attribute2)]][attribute1_restrictions]
        if (is.null(x$condition[[1]])) {
          cat(paste0("\nConditional on ", attribute2, " in {", paste(conditions_cacp1, collapse = ", "), "}:\n"))
          print.default(format(x$estimates[x$indices$cacp_byattr[[attr_nb]]$indices], digits = digits), print.gap = 2L, quote = FALSE)
        }
        if (!is.null(x$condition[[1]])) {
          if (attribute1 != x$condition[[1]]) {
            cat(paste0("\nConditional on ", attribute2, " in {", paste(conditions_cacp1, collapse = ", "), "}:\n"))
            print.default(format(x$estimates[x$indices$cacp_byattr[[attr_nb]]$indices], digits = digits), print.gap = 2L, quote = FALSE)
          }
        }
      }
    } else if (cacp == "comparable_pairs") {
      cat("\nConditional on comparable pairs:\n")
      print.default(format(x$estimates_alt[x$indices$cacp], digits = digits), print.gap = 2L, quote = FALSE)
    }
  } else {
  cat("No coefficients\n")
  cat("\n")
  invisible(x)
  }
}





### SUMMARY ###########################################################


summary.conjacp <- function(object,
                            cacp = "conditional") {
  z       <- object
  se      <- sqrt(diag(z$vcov))
  est     <- z$estimates
  se_alt  <- sqrt(diag(z$vcov_alt))
  est_alt <- z$estimates_alt
  ans     <- z[c("subgroups_parameters", "n_profiles", "n_tasks", "n_respondents", "subset",
                 "estimand", "attr_list", "attr_restricted", "levels_restricted",
                 "attr_levels", "condition", "indices")]
  if (object$estimand %in% c("acp", "p")) {
    ans$coefficients     <- cbind(est, se)
    ans$coefficients_alt <- cbind(est_alt, se_alt)
    dimnames(ans$coefficients)     <- list(names(z$estimates), c("Estimate", "Std. Error"))
    dimnames(ans$coefficients_alt) <- list(names(z$estimates_alt), c("Estimate", "Std. Error"))
  } else {
    ans$coefficients     <- cbind(est, se, est/se, 2*pnorm(abs(est/se), lower.tail = FALSE))
    ans$coefficients_alt <- cbind(est_alt, se_alt, est_alt/se_alt, 2*pnorm(abs(est_alt/se_alt), lower.tail = FALSE))
    dimnames(ans$coefficients)     <- list(names(z$estimates), c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
    dimnames(ans$coefficients_alt) <- list(names(z$estimates_alt), c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
  }
  ans$cacp   <- cacp
  class(ans) <- "summary.conjacp"
  ans
}

print.summary.conjacp <- function (x,
                                   digits = max(3L, getOption("digits") - 3L)
)  {
  cat("\nNumber of profiles:", paste0(x$n_profiles))
  cat("\nNumber of tasks:", paste0(x$n_tasks))
  cat("\nNumber of respondents:", paste0(x$n_respondents))
  if (length(x$indices$cacp_byattr)) {
    cat("\nRestrictions:")
    if (length(x$levels_restricted) == 2) {
      attribute1              <- x$levels_restricted[[1]]$attribute1
      attribute2              <- x$levels_restricted[[1]]$attribute2
      attribute1_restrictions <- as.numeric(unlist(x$levels_restricted[[1]]$restrictions_attribute1))
      attribute1_restrictions <- x$attr_levels[[which(x$attr_list %in% attribute1)]][attribute1_restrictions]
      attribute2_restrictions <- as.numeric(unlist(x$levels_restricted[[1]]$restrictions_attribute2))
      attribute2_restrictions <- x$attr_levels[[which(x$attr_list %in% attribute2)]][attribute2_restrictions]
      cat("\n ", paste0("   [", attribute1, "] ",
                        paste(attribute1_restrictions, collapse = ", "),
                        "\n   & [", attribute2, "] ",
                        paste(attribute2_restrictions, collapse = ", ")))
    } else {
      for (i in 1:(length(x$levels_restricted)/2)) {
        k <- 2 * i - 1
        attribute1              <- x$levels_restricted[[k]]$attribute1
        attribute2              <- x$levels_restricted[[k]]$attribute2
        attribute1_restrictions <- as.numeric(unlist(x$levels_restricted[[k]]$restrictions_attribute1))
        attribute1_restrictions <- x$attr_levels[[which(x$attr_list %in% attribute1)]][attribute1_restrictions]
        attribute2_restrictions <- as.numeric(unlist(x$levels_restricted[[k]]$restrictions_attribute2))
        attribute2_restrictions <- x$attr_levels[[which(x$attr_list %in% attribute2)]][attribute2_restrictions]
        cat("\n ", paste0(i, ".   [", attribute1, "] ",
                          paste(attribute1_restrictions, collapse = ", "),
                          "\n     & [", attribute2, "] ",
                          paste(attribute2_restrictions, collapse = ", ")))
      }
    }
  }
  if (!is.null(x$subset)) {
    cat("\nSubset:", paste0(x$subset[[1]], " == ", x$subset[[2]]))
  }
  if (!is.null(x$condition)) {
    cat("\nCondition:", paste0(x$condition[[1]], " in {", paste(x$condition[[2]], collapse = ", "), " }"))
  }
  if (x$estimand == "dacp") {
    cat("\nSubgroups:", paste0(x$subgroups_parameters$var))
  }
  cat("\n\n")
  if (length(x$indices$acp)) {
    if (x$estimand == "acp") {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nACPs:\n")
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nACPs: (unrestricted attributes)\n")
      }
    } else if (x$estimand == "dacp") {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nDACPs,",
            paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":\n"))
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nDACPs,",
            paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":  (unrestricted attributes)\n"))
      }
    } else {
      if (length(x$indices$cacp_byattr) == 0) {
        cat("\nDirect Pairwise Preferences:\n")
      } else if (length(x$indices$cacp_byattr)) {
        cat("\nDirect Pairwise Preferences: (unrestricted attributes)\n")
      }
    }
    coefs_acp <- x$coefficients[x$indices$acp,]
    printCoefmat(coefs_acp, digits, na.print = "NA")
  }
  if (length(x$indices$cacp_byattr)) {
    if (x$estimand == "acp") {
      cat("\nCACPs:\n")
    } else if (x$estimand == "dacp") {
      cat("\nDCACPs,",
          paste0(x$subgroups_parameters$ref, " - ", x$subgroups_parameters$oth, ":\n"))
    } else {
      cat("\nDirect Pairwise Preferences:\n")
    }
    if (x$cacp == "conditional") {
      for (attr_nb in 1:length(x$indices$cacp_byattr)) {
        attribute1              <- x$indices$cacp_byattr[[attr_nb]][[1]]
        attribute2              <- x$levels_restricted[[which(x$attr_restricted %in% attribute1)]]$attribute2
        attribute1_restrictions <- as.numeric(x$levels_restricted[[which(x$attr_restricted %in% attribute1)]]$restrictions_attribute2)
        conditions_cacp1        <- x$attr_levels[[which(x$attr_list %in% attribute2)]][-attribute1_restrictions]
        conditions_cacp2        <- x$attr_levels[[which(x$attr_list %in% attribute2)]][attribute1_restrictions]
        if (is.null(x$condition[[1]])) {
          cat(paste0("\nConditional on ", attribute2, " in {", paste(conditions_cacp1, collapse = ", "), "}:\n"))
          coefs_cacp1 <- x$coefficients[x$indices$cacp_byattr[[attr_nb]]$indices,]
          printCoefmat(coefs_cacp1, digits, na.print = "NA")
        }
        if (!is.null(x$condition[[1]])) {
          if (attribute1 != x$condition[[1]]) {
            cat(paste0("\nConditional on ", attribute2, " in {", paste(conditions_cacp1, collapse = ", "), "}:\n"))
            coefs_cacp1 <- x$coefficients[x$indices$cacp_byattr[[attr_nb]]$indices,]
            printCoefmat(coefs_cacp1, digits, na.print = "NA")
          }
        }
      }
    } else if (x$cacp == "comparable_pairs") {
      cat("\nConditional on \"comparable pairs\":\n")
      coefs_cacp2 <- x$coefficients_alt[x$indices$cacp,]
      printCoefmat(coefs_cacp2, digits, na.print = "NA")
    }
  } else {
  cat("\n")
  cat("\n")
  invisible(x)
  }
}





### CONJACP.DIFF ###########################################################


conjacp.diff <- function(conjacp.object,
                         acp1,
                         acp2,
                         cacp = "conditional",
                         digits = max(3L, getOption("digits") - 3L),
                         signif.stars = getOption("show.signif.stars")
) {
  if (cacp == "conditional") {
    estimates <- conjacp.object$estimates
    vcov      <- conjacp.object$vcov
  } else if (cacp == "comparable_pairs") {
    estimates <- conjacp.object$estimates_alt
    vcov      <- conjacp.object$vcov_alt
  }
  index1    <- which(names(estimates) %in% acp1)
  index2    <- which(names(estimates) %in% acp2)
  diff      <- estimates[index1] - estimates[index2]
  se        <- sqrt(vcov[index1,index1] + vcov[index2,index2] - 2 * vcov[index1,index2])
  zval      <- diff/se
  out       <- cbind(diff, se, zval, 2*pnorm(abs(zval), lower.tail = FALSE))
  dimnames(out) <-
    list(paste0(names(estimates[index1]), "-", names(estimates[index2])),
         c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
  printCoefmat(out, digits, signif.stars = signif.stars,
               na.print = "NA")
}





### CONJACP.VAR ###########################################################


conjacp.var <- function(conjacp.object,
                        cacp = "conditional",
                        alpha = .05,
                        nsimul = 1000) {
  
  ## Packages
  
  require(MASS)
  require(stringr)
  "%nin%" <- Negate("%in%")
  
  
  ## Transition from conjacp.data
  
  # Get data from conjacp.object
  if (cacp == "conditional") {
    estimates       <- conjacp.object[[1]]
    vcov            <- conjacp.object[[2]]
  } else if (cacp == "all_comparable") {
    estimates       <- conjacp.object[[3]]
    vcov            <- conjacp.object[[4]]
  }
  estimand        <- conjacp.object[[5]]
  attr_list       <- conjacp.object[[11]]
  attr_restricted <- conjacp.object[[12]]
  condition       <- conjacp.object[[15]]
  indices         <- conjacp.object[[16]]
  
  # Sanity check
  if (estimand != "acp") {
    stop("Within-attribute variability estimates only available for ACPs.")
  }
  
  
  ## Simulate posterior distribution
  
  # Create empty matrices to store simulations
  range_estimates_s       <- matrix(NA,
                                    nrow = nsimul,
                                    ncol = length(attr_list[which(attr_list %nin% condition[[1]])]))
  variability_estimates_s <- matrix(NA,
                                    nrow = nsimul,
                                    ncol = length(attr_list[which(attr_list %nin% condition[[1]])]))
  
  # Draw ACP estimates from their posterior distributions
  estimates_posdef <- estimates[duplicated(str_extract(names(estimates), ".+?(?=\\.)")) == TRUE]
  vcov_posdef <- vcov[duplicated(str_extract(names(estimates), ".+?(?=\\.)")) == TRUE,
                      duplicated(str_extract(names(estimates), ".+?(?=\\.)")) == TRUE]
  estimates_s <- mvrnorm(nsimul, estimates_posdef, vcov_posdef)
  
  # Calculate range and variability distributions
  for (attr_temp in attr_list[which(attr_list %nin% condition[[1]])]) {
    
    attr_nb                         <- which(attr_list[which(attr_list %nin% condition[[1]])] %in% attr_temp)
    index                           <- which(str_extract(names(estimates_posdef), ".+?(?=\\.)") == attr_temp)
    
    if (length(index) == 1) {
      range_estimates_s[,attr_nb]       <- 2 * abs(estimates_s[,index])
      variability_estimates_s[,attr_nb] <- abs(estimates_s[,index])
    } else if (length(index) > 1) {
      estimates_s_temp                  <- cbind(estimates_s[,index], as.matrix(-apply(estimates_s[,index], 1, sum)))
      range_estimates_s[,attr_nb]       <- apply(estimates_s_temp, 1, function(x) max(x) - min(x))
      variability_estimates_s[,attr_nb] <- apply(estimates_s_temp, 1, function(x) mean(abs(x)))
    }
    
  }
  
  
  ## Calculate QOIs
  
  labels <- attr_list[which(attr_list %nin% condition[[1]])]
  
  # Range
  range_estimates        <- apply(range_estimates_s, 2, mean)
  names(range_estimates) <- labels
  range_vcov             <- cov(range_estimates_s)
  dimnames(range_vcov)   <- list(labels, labels)
  range_lower            <- apply(range_estimates_s, 2, function(x) quantile(x, .025))
  names(range_lower)     <- labels
  range_upper            <- apply(range_estimates_s, 2, function(x) quantile(x, .975))
  names(range_upper)     <- labels
  
  # Variability
  variability_estimates        <- apply(variability_estimates_s, 2, mean)
  names(variability_estimates) <- labels
  variability_vcov             <- cov(variability_estimates_s)
  dimnames(variability_vcov)   <- list(labels, labels)
  variability_lower            <- apply(variability_estimates_s, 2, function(x) quantile(x, .025))
  names(variability_lower)     <- labels
  variability_upper            <- apply(variability_estimates_s, 2, function(x) quantile(x, .975))
  names(variability_upper)     <- labels
  
  
  ## Out
  return(list(range_estimates       = range_estimates,
              range_vcov            = range_vcov,
              range_lower           = range_lower,
              range_upper           = range_upper,
              variability_estimates = variability_estimates,
              variability_vcov      = variability_vcov,
              variability_lower     = variability_lower,
              variability_upper     = variability_upper))
  
  
}
  
